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Abstract 

We investigate the stability properties of two different classes of metabolic cycles using a combination 
of analytical and computational methods. Using principles from structural kinetic modeling (SKM), 
we show that the stability of metabolic networks with certain structural regularities can be studied 
using a combination of analytical and computational techniques. We then apply these techniques 
to a class of single input, single output metabolic cycles, and find that the cycles are stable under 
all conditions tested. Next, we extend our analysis to a small autocatalytic cycle, and determine 
parameter regimes within which the cycle is very likely to be stable. We demonstrate that analytical 
methods can be used to understand the relationship between kinetic parameters and stability, and that 
results from these analytical methods can be confirmed with computational experiments. In addition, 
our results suggest that elevated metabolite concentrations and certain crucial saturation parameters 
can strongly affect the stability of the entire metabolic cycle. We discuss our results in light of the 
possibility that evolutionary forces may select for metabolic network topologies with a high intrinsic 
probability of being stable. Furthermore, our conclusions support the hypothesis that certain types 
of metabolic cycles may have played a role in the development of primitive metabolism despite the 
absence of regulatory mechanisms. 



1. Introduction 

Cycles are at the heart of the metabolic networks of organisms spanning the entire tree of life 
(USUI- For example, the tricarboxylic acid (TCA) cycle sits at the core of energy production for many 
species and additionally plays the role of regenerating essential cellular nutrients and components. A 
great deal of research has been devoted to understanding the stability properties of such cycles, which 
have often been represented as chemical reaction networks (CRNs). Many of these studies have drawn 
conclusions about the steady state of a CRN without regard to details about the rate constants or 
chemical concentrations of the particular CRN (for instance, see [3, [H, [|| ) . Instead, these approaches 
rely on the topology of the network and commonly assume mass-action kinetics laws to limit (and 
sometimes determine) the possible behaviors the system may exhibit at steady state. Can one hope 
to extend these results to general, potentially nonlinear, rate laws? 

In this work, we use a generalized modeling framework known as structural kinetic modeling (SKM) 
to extend previous results on two classes of metabolic cycles. The first class has a special but relatively 
flexible structure: it contains a single input, a single output, and unlimited length. This relatively 
organized structure makes it amenable to analytical methods of study. For cycles belonging to this first 
class, we prove that the cycle can only lose dynamical stability in an oscillatory manner, regardless of 
absolute metabolite concentrations, flux values, kinetic rate constants, and form of monotonic kinetic 
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rate law. We further investigate the cycle computationally, and find that it remains completely stable 
under all conditions tested, leading us to make a conjecture about its stability under all conditions. 

We then apply the same SKM approach to a second class of cycles, autocatalytic cycles, and 
demonstrate that these are not stable under all combinations of kinetic parameters and metabolic 
conditions. Instead, we use analytical calculations to derive relationships between key parameters 
which can enhance or attenuate the stability properties of the cycle. These relationships are then 
computationally tested and verified. 

Nearly all prior published work using SKM @,@] employs computational experiments and statistical 
methods in order to draw conclusions about the system of interest. In this work, we show that metabolic 
systems with relatively well-organized structures are not only analytically tractable but furthermore 
that their analysis can actually lead to interesting and novel conclusions. Finally, we discuss the 
implications of our results on the evolution of primitive metabolic networks. The SKM technique 
used throughout the paper highlights essential parameters which directly impact the stability of the 
cycle under study. Particular combinations of these parameters may have conferred an evolutionary 
advantage to their respective metabolic networks by making them robust to small environmental 
perturbations. 



2. Background 

2.1. Structural Kinetic Modeling (SKM) 

The methods developed in t his p aper are based on the SKM framework^. SKM is a specific appli- 
cation of generalized modeling [10] in which normalized parameters replace conventional parameters 
(such as V max or Km in the modeling of metabolic networks). The normalized parameters have a 
direct connection to the original kinetic parameters, but are much easier to work with. As will be 
shown below, these parameters usually have well-defined and limited ranges (e.g. [0,1]), and sampling 
them across this range effectively samples all possible values of the original kinetic parameters. 

The goal of SKM is to capture the local stability properties of a biochemical system. In this sense, 



it bridges the gap between genome-scale steady state modeling [ll| and explicit kinetic modeling of a 
metabolic process. To study stability, one usually determines the Jacobian J of the system of interest 
and evaluates it at the system's steady states. Assuming knowledge of the form of kinetic rate laws 
for a CRN, it is quite possible to write down the corresponding J. However, in many cases J will be 
unnecessarily complicated and quite difficult to work with. By performing a change of variables, SKM 
actually simplifies the mathematical form of each entry in J. The new entries are in almost all cases 
easier to work with. The specific method by which the change of variables occurs is briefly described 



below and more completely in | Appendix A Much of this material is based on Q and we refer readers 



to that reference and its supplementary materials for more information. 

If we let S be the m-dimensional vector of metabolite concentrations, N be the mxr stoichiometric 
matrix, and v be the r-dimensional vector of reaction rates, then we can describe the dynamics of the 
system with the equation 

H = N»(S,k) (2.1) 

where u(S,k) denotes that the reaction rates are dependent on both metabolite concentrations S 
and kinetic parameters (such as Michaelis-Menten constants) k. If we assume that a non-negative 
steady state S° exists, then we can redefine our system using the definitions 

,-Ma.^^ Mi) -^ (22) 
x 'i — go i-'mj — JV y go >/M x ) — ^j. (S°) 

where i = l...m and j = l...r, x is a vector of metabolite concentrations normalized with respect 
to their steady state concentrations and fi represents flux normalized with respect to steady state flux 
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values. The matrix A represents the stoichiometric matrix normalized with respect to steady state 
fluxes and steady state metabolite concentrations. 

The Jacobian, evaluated at x = 1 (which, because of the way x is defined, is actually the equilib- 
rium of the system) , can be written as 



J x = A0£ (2.3) 

Note that the equations above were derived without regard to the actual form of the kinetic 
equations that determined the ODE system. The matrix 9£ contains elements which represent the 
degree of saturation of normalized flux fij with respect to normalized substrate concentration Xj. 
In terms of derivatives, each element of 9 represents the degree of change in a flux as a particular 
metabolite is incrementally increased. This is analogous to the concept of elasticity in metabolic 
control analysis [12j . 

What does the 9 matrix look like? Its columns correspond to each metabolite, and its rows to each 
flux. A non-zero element 9* in the matrix represents the effect a small change in a metabolite j has on 
flux i. In the case of Michaelis-Menten kinetics, this element in the matrix may take values ranging 
from [0,1]. In the case of standard competitive inhibition (e.g. allosteric inhibition by a product), the 
element takes values in [-1,0]. 

To usefully illustrate the meaning of a single 9 element, consider an equation following Michaelis- 
Menten kinetics shown in (12.41) . First, we write n{x) 1 which we recall is the flux normalized by the 
flux at the steady state. Here, Sq is the concentration of the substrate S at steady state. Then, we 
manipulate the equation by substituting xSq for S, where x is the normalized steady state concentration 
of substrate S. The result is shown in ()2.5|) . 

k 2 E S T s Q 

Finally, we take a derivative with respect to x and evaluate it at x — 1 to obtain 9 in (12.61) . Notice 
that 9 can only take values between and 1 for any positive value of Sq . 
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(2.6) 



The power of SKM comes from the parametrization illustrated above. Each element in 9 has a 
precise correspondence to some combination of kinetic parameters in the original model. However, it 
is far more tractable to study a system using 9 parameters rather than the original kinetic parameters. 
To see this more clearly, recall that in most cases, biochemical kinetic constants are poorly estimated. 
In order to build precise ODE models for biochemical systems, it is usually necessary to actually choose 
values for these constants. While the chosen values may be estimated from experimental measurements, 
hidden dependencies in the model may actually result in non-obvious correlations between parameters 
that can strongly affect the output of the model. 

On the other hand, if one is only interested in the stability of a characterized steady-state, it may 
not be necessary to actually have knowledge of kinetic parameters. Experimental measurements can 
provide data on absolute metabolite concentrations and flux values, and the stoichiometry is often 
known a priori. Then, one can parameterize the system as shown above and sample many possible 
combinations of 9 parameters. For each unique set of 9 parameters, the stability of the Jacobian is 
determined and recorded. Analysis of the results can lead to several important conclusions such as 
which 9 parameters have the strongest correlation with stability of the entire metabolic system. Thus, 
the benefit of employing SKM over other techniques is that it provides the means to analyze and make 
sense of a large number of possible cases of a metabolic network, rather than a single instance. 
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At first glance, it is unclear whether this will make things easier or harder; values for fluxes and 
concentrations are almost never known. However, our work to follow shows that the absence of this 
knowledge does not make the problem of stability intractable. Instead, we operate with the following 
idea: we assume that the steady state flux vector v and the metabolite concentration vector X are 
variables, and we do not specify them. Instead, we simply determine the stability of the system in 
terms of these and other (9) variables. Then, we can find trends in stability as the flux and metabolite 
concentrations change: it may be that as a concentration goes up, the stability of the system tends to 
fall. Isolating these trends is the motivation of our analysis. 



2.2. Previous results on stability of metabolic cycles 

Our work is preceded by a great deal of work aimed at understanding the stability of metabolic 
networks [1 M, 0, El E3] . One flavor of such work is referred to as chemical reaction network theory 
(CRNT). CRNT is an elegant and powerful framework for understanding how the structure of chemical 
networks dictates their equilibrium and stability properties. By applying CRNT to the cycles studied 
in this paper, one can (in some cases) prove the stability of their steady states. This result can be 
easily obtained, for example, by feeding the stoichiometric matrices of the corresponding networks 
into the ERNEST MATLAB Toolbox [13] ■ Notably, these results cannot be easily extended beyond 
the case of mass-action kinetics (where the rate of reaction is proportional to the concentration of its 
reactants). The SKM method presented in this article extends the stability results for these cycles to 
general monotonic kinetics. 

Dynamical systems approaches have also been used in understanding the stability properties of 



autocatalytic cycles (the second type of cycle presented in the paper) . Previous research by King 14 [ 
(focusing on dynamics) and by Eigen |3( (focusing on steady-state stability properties) have studied 
autocatalytic cycles assuming first-order (non-saturating) mass-action kinetic rate laws. Extending 
such work, our analysis takes into account all forms of saturating rate laws up to order one, and 
derives simple expressions for the stability of metabolic cycles given adequate information about rate 
constants and steady-state concentrations. To put this in context, the prior work on autocatalytic 
cycles represents a single point (at which all 9 parameters equal 1) in the large space of 9 parameters 
within which the actual metabolic network resides. Our approach searches over this entire space, and 
identifies trends in stability across it. 



3. A Simple Example 

To illustrate the main concepts behind analytically characterizing the stability of a metabolic cycle 
using SKM, we will start with a simple example. Consider the cycle shown in Figure [T] To analyze 
the conditions under which it is stable, we will write the A and 9 matrices, multiply them together to 
obtain the Jacobian </, and then find the eigenvalues of J by writing out its characteristic equation. 

The cycle contains 3 metabolites and 6 reactions. We impose that the reactions can only proceed 
in the forward direction. Furthermore, the 6 reactions are constrained by mass balance, and we 
note that two linearly independent rate vectors will characterize the fluxes of all 6 reactions in the 
network. Assume that the steady states metabolite concentrations (A , B°, C°, O®, 0%) °f the system 
are (1,1,1,1,1). Furthermore, assume that a flux of magnitude aF enters the cycle through reaction 
Vi, and (1 — a)F flux returns through reaction U5 (where < a < 1). 

Finally, we impose that there is no complicated activation or inhibition present in the system. This 
means that vi uses only metabolite A and cofactor 0\ as its substrates, V3 uses only metabolite B, V4 
and U5 use only metabolite C, vq uses only 2 , and v\ is constant and uses none of the metabolites 
as substrates. Notice that this constrains all 9 parameters to be greater than or equal to zero. Con- 
straining 9 to be positive ensures that any kinetic rate laws we consider in our analysis are monotonic, 
i.e. that an increase in the concentration of the substrate of a reaction will never decrease the rate of 
a reaction. This will be essential in the analysis to follow. Now, we can write our system in terms of 
SKM variables: 
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Figure 1: A simple metabolic cycle 
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(3.2) 



Because the cofactors 0\ and Oi are conserved, we can remove them from the two matrices above 
(for details and justification, see the Supplementary Material of [9|). Doing this, we omit the last row 
of A, and 9 becomes 



(3.1) 







(3.3) 



Then, to obtain the Jacobian we compute the product A6. In order to determine the stability of 
the system, we find the eigenvalues of the Jacobian and determine whether any of them are positive. 
We hnd the eigenvalues by imposing | J — XI \ = 0, i.e. by setting to zero the following determinant: 



-A - F9i 

F6\ — A — FO2 

F6» 2 

-6» 2 

This corresponds to the equation 
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{-o0 3 - (1 - a)9 4 - A)(0i - A) ((-9 5 - 9 6 - X)(-9 2 - A) - 9 2 9 5 ) + 

(l-a)9 1 e 2 e 4 (-X-9 6 ) = 



(3.5) 
(3.6) 



The stability properties will be determined by the signs of the real parts of the eigenvalues obtained 
from solving this equation. For all eigenvalues with only real components, it is possible to determine 
the signs without explicitly solving the equation, by using Descartes' Rule of Signs. If we arrange the 
characteristic polynomial in order of decreasing variable exponents, then the number of positive real 
roots is equal to the number of signs changes in the coefficients of the ordered polynomial, or less than 
the number of size changes by a multiple of 2. In our case, there are no sign changes (all coefficients 
are greater than zero, notice that any negative terms cancel), so there are no positive roots. Note 
that it is not possible to have a zero eigenvalue either, because there is a nonzero coefficient for the A 
term. Then, we must consider two alternative cases. In one, the Jacobian has only real eigenvalues, 
in which case they are all negative and the system is stable. In the other, the Jacobian may contain 
pairs of complex conjugate eigenvalues whose real part will determine the stability of the system. To 
investigate systems which may display such complex eigenvalues, we sampled many instances of the 
Jacobian computationally and computed the corresponding eigenvalues. Any eigenvalues with real 
part greater than zero would indicate instability. For this particular cycle, we randomly sampled each 
9 in the range (0,1] for a value of a in the open set (0,1). Then, the eigenvalues of the Jacobian were 
calculated. This process was repeated 4 x 10 6 times and in no cases were any unstable Jacobians 
identified. 



4. A More General Treatment 

In this section we generalize our previous results to a metabolic cycle of arbitrary length, relying 
on the following assumptions: 

1. There exists at least a single steady state for the system. It is trivial to show that such a 
steady state exists because rate constants (such as V max ) are not related to any 9 parameters 
which determine the stability, and can be chosen independently and arbitrarily in order to satisfy 
the steady-state condition. Furthermore, results from [4| and computational analysis using the 



ERNEST MATLAB package [17( demonstrate that the cycle under study is strongly sign deter- 
mined (SSD) and injective. It therefore admits only a single unique equilibrium for a given set 
of kinetic parameters (e.g. rate constants, saturation constants, etc.). It is worthwhile to note 
that if multiple equilibria did exist, the results presented in this section would still be valid. 
The input and output reactions from the cycle are separated by one metabolite. In other words, if 
the output reaction of the cycle leaves through metabolite Xi in the cycle, then the input reaction 
to the cycle enters through metabolite Xi+x. In analogy to the simple example above, the output 
reaction leaves through metabolite C, and the input reaction enters in the next metabolite in 
the chain, A. This assumption is made to simplify the analysis, but computational experiments 
indicate it does not appear to be necessary. However, this has not been explicitly proven. 
The j th reaction (j = 2, 3, N + l) uses only the (j — l) th metabolite as a substrate. There is no 
feedback inhibition or activation of any kind. This means that all 9 parameters are constrained 
to be greater than zero. A negative 9 would correspond to a metabolite having an inhibitory 
effect on the rate of reaction. 

Cofactors such as ATP and NADH are often used as the energetic driving forces in metabolic 
reactions in order to make these reactions thermodynamically feasible. We assume that a single 
cofactor pair (0\ and 2 ) is involved in the metabolic cycle, and that it is not involved in the 
reaction creating the first or last metabolite in the cycle. The reformation of cofactors is assumed 
to be driven by some energy input into the system. 
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Figure 2: A complete metabolic cycle with cofactors 



Our goal is to study the stability of the cycle shown in figure [U which generalizes the cycle of 
Figure [TJ 

The generalized form of A and 8, analogous to the ones from (|3.1[) and (|3.3p . is shown in Appendix 
|Appcndix C Note that we make no further assumptions about the steady states of the system: 
X° — (Xi,X2, ■ ■ ■ , A^) and the flux through the system remains equal to F. With N metabolites 
and N + 3 reactions, A is N x (TV + 3) and 9 is (N + 3) x N. The complete analysis is shown in 
|Appcndix C[ and the results are summarized below. 

The characteristic equation for the eigenvalues of the system reduces to 
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The terms on the bottom two lines will cancel with terms from the top two lines, resulting in a 
polynomial equation of order N + 1 containing all nonzero coefficients with exactly the same sign. 
Furthermore, a constant A term will remain, indicating that there is no zero eigenvalue. Therefore, 
all eigenvalues with only real part are negative, and the system can only lose stability in an oscillatory 
manner. 

We proceeded with computational experiments in which SKM parameters were randomly sampled 
to search for cases of systems whose eigenvalues lie in the positive real part of the complex plane. In 
these experiments, we generated random samples of 9 parameters (in the range (0,10] to account for 
high Hill coefficients), a in (0,1), metabolite concentrations in (0,100), and calculated the eigenvalues 
of the Jacobian. Approximately 4 x 10 6 different instances were simulated for each particular cycle 
length from 3 to 8, and for each possible cofactor location. In none of the sampled instances was the 
real part of any eigenvalue of the Jacobian greater than zero, indicating that all samples were stable. 

In a separate computational experiment, we tested the dependency of our stability results on the 
structure of the cycle. An example is shown in the lower part Figure [3] in which the first metabolite 
activates the reaction creating the third metabolite (similar to a feedforward loop). This introduces 
an instability into the system which does not appear in our prototypical system. 

Thus, both computational experiments and analytical results support the conjecture that under 
general kinetic rate laws, cycles of this form are stable under all conditions. We believe that a proof of 
this conjecture may follow from methods similar to those used to establish stability results for linear 
chains with positive feedback, as presented in [1 61 ] - 

5. Autocatalytic Cycles 

The next question we want to address is whether the methods and results we derived for the 
cycles of metabolic reactions seen in the previous sections can be extended to cycles of fundamentally 
different topology, namely autocatalytic cycles. Autocatalytic cycles consist of a set of metabolites 
which reproduce themselves with each turn of the cycle. Their simplified stoichiometry often takes the 
form A + X — > 2A, where a substrate X transforms into a product A, through catalysis of A itself 
[II EH- The formose reaction and the reductive citric acid cycle are two cases of autocatalytic cycles; 
both have been hypothesized to play roles in early metabolism, and significant work has been devoted 
to understanding the conditions under which they and other autocatalytic cycles may have naturally 
appeared and proliferated @, . 

Shown in Figure [4j our simple example contains 3 metabolites and 5 reactions. It is distinct from 
the cycles shown in prior sections because both the first and last metabolite (A and C) are required 
for reaction V2 to occur. Furthermore, we assume that reactions exiting from the cycle transport 
metabolites to a separate compartment, where they are not accessible to the reactions within the 
cycle. Notably, we will show that this system does not contain a universally stable steady state. 
Instead, the steady state remains stable under certain parameter regimes, to be defined in the analysis 
below. Although this example is relatively simple, the results can be extended to more complicated 
examples using a method similar to the one described in Appendix C. 

The same assumptions are made as before, and all reactions are assumed to proceed in only the 
forward direction. There are 2 independent flux vectors which span the null space of the stoichiometric 
matrix for this system. Therefore, we assign the flux vector V — (vi,V2,Va,V4,,Vs) = (aF, F, F(x — 
7), jF, F(x — 7 — 1)), where jF is the outflow from the cycle through the reaction V4 and F(x — 7 — 1) 
is the outflow from the cycle through the reaction U5. Furthermore, it is noted that for flux balance 
considerations, x > 7. Since all reactions are assumed to go in the forward direction, x > 7 + 1. 

Following the same analysis as before, we can write the matrices A and 0: 
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Figure 3: | (a) | Histogram of values of the maximal eigenvalue for a case of the system shown in Figure [2] The specific 
system contains 6 metabolites, and the reaction vi is coupled to a cofactor. All eigenvalues are stable, indicating the 
system is stable. |(b)| Histogram of values of the maximal eigenvalue for a system in which one of the assumptions is 
broken. In this case, metabolite X\ activates reaction 113. Notice that not all maximal eigenvalues are negative, indicating 
that the system may be unstable for certain combinations of parameters. This suggests that global stability strongly 
relies on the topology of the network, and that additional activating or inhibitory reactions may prevent a steady state 
from being stable. 
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Then, the Jacobian J is the product of A and 8: 
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Figure 4: A simple auto-catalytic cycle 



By imposing | J — \I\ = 0, we get the characteristic equation 
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Expanding this equation, the coefficients of the third-order polynomial in A are: 
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To analyze the signs of the roots of this polynomial, we use the Routh-Hurwitz criterion. According 
to this criterion, a third-order polynomial in z with coefficients, a, 6, c of the form z 3 + az 2 + bz + c = 
has roots lying in the negative real half of the complex plane if the conditions a > 0,c> 0, ab — c > 
are satisfied. For our polynomial, it is obvious that both a and c are always greater than zero as long as 
x > 7+ 1, so we are left to find conditions where the last inequality is satisfied. Note that if x < 7 + 1, 
then c < 0, there is a positive root to the characteristic equation, and the system is unstable. 

Rather than deriving a precise condition for stability in terms of any single variable, we can use 
simple observations regarding the particular forms of the coefficients to infer trends in stability. In 
particular, notice that all negative terms in — c are precisely canceled out by corresponding terms in 
ab. Then, the only term which may cause ab — c to be less than zero is (^-t^M 1 -^) ^ ^here are 
several ways to minimize the negative value of this term, such as: decreasing x (while ensuring that 
x > 7 + I), increasing 7 (while ensuring that x > 7 + 1), decreasing 02 and/or 04, increasing 6\ and/or 
03 and/or 5 , decreasing A , and increasing B or Co- 

We can computationally test these predictions by following a procedure similar to Q and ran- 
domly sampling the parameters that define the autocatalytic cycle. Specifically, we can calculate the 
eigenvalues of the Jacobian matrix as a single parameter of interest (such as a parameter) is varied, 
while all other parameters and metabolite concentrations are sampled randomly. By doing so, we 
sample over many possible instances of the metabolic network, and can identify correlations between 
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individual parameters and stability. For each combination of dependent variable (e.g.: x, 61,62, Aq) 
and 7 shown in the figures below, 5 x 10 5 instances of the Jacobian were sampled and characterized. 
The value of the maximal eigenvalue for each instance indicated whether the network is stable or not. 

From the terms of the A coefficient derived above, one can see that as long as x > 7 + 1, it is 
ideal for 7 to be as large as possible and x to be as small as possible. This minimizes the value of 
(£-7)^2^(1-3;) ^ j n pjg ure 5( a )| stability of the system is plotted against values of the stoichiometric 
coefficient x of metabolite B in reaction V2- Each line indicates a different value for the outflow 7. As 
x increases, the stability of the system falls. Additionally, as the value of 7 increases, the stability of 
the system increases, while the general shape of the plot remains identical. Note that if a; < 7 + 1, then 
the system is universally unstable, in agreement with the condition derived before. It is particularly 
interesting that the system's optimally stable operating point is located at the boundary of a region 
where it is entirely unstable due to dilution effects. 





(a) 



(b) 





(c) (d) 

Figure 5: Each point in the figure is determined by sampling the Jacobian 5 X 10 5 times, [(a)] Proportion of stable steady 
states versus stoichiometric coefficient x, for different values of 7, |(b)| versus 6\, for different values of 7, |(c)| versus 82, 
for different values of 7, |(d)|versus Aq, for different values of 7 



Similarly, Figure 5(b) confirms that as 6\ increases, the stability of the system rises. In Figure 5(c) 
it is shown that as 62 increases, the stability of the system falls. Furthermore, the relationship between 
stability and 7 is retained in all these plots. 

The combination of analytical and computational results presented so far in this section indicates 
that the tuning of a single parameter can greatly increase the likelihood of a metabolic system being 
stable. One may wonder which parameters, among those that can affect stability, are the ones that 
a cellular system would be most likely to adjust. An intriguing possibility supported by the current 
findings is that a cell could approach a stable steady state by appropriately tuning metabolite con- 
centrations. A sufficiently high steady state concentration of B or C would nearly guarantee that the 
simple autocatalytic cycle would be stable. In the future it may be interesting to verify experimentally 
whether certain metabolic pathways - and autocatalytic cycles in particular - may display steady state 
metabolite concentrations that fall into the range of high stability. 
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6. Discussion 



We have used structural kinetic modeling (SKM) to assess the dynamical stability of certain classes 
of metabolic cycles. In prior work, it was demonstrated that SKM can be used as a powerful computa- 
tional tool for identifying stabilizing sites in metabolic networks by sampling large ranges of parameters 
Q ■ The main thrust of our work was to show that SKM can also be used analytically to characterize 
stability. Drawing from basic concepts in nonlinear dynamics and linear algebra, we demonstrated that 
it is possible to make a strong assertion about the stability of a class of non-autocatalytic metabolic 
cycles, and to identify trends in stability as a function of model parameters for a type of autocatalytic 
cycles. 

For the non-autocatalytic cycles of the type illustrated in Figure [3J we have proven that the topol- 
ogy of the networks ensures that they can only lose stability in an oscillatory manner, regardless of 
metabolite concentrations, flux magnitude, or form of kinetic equations. The generality of the result 
with respect to kinetic laws is particularly noteworthy. Computational experiments support the con- 
jecture that these cycles are completely stable and do not produce Jacobians with positive complex 
eigenvalues. The stability of non-autocatalytic cycles similar to the one studied in this section could 
be easily proven using the Deficiency One Theorem of Chemical Reaction Network Theory (CRNT), 
However the CRNT result, in this case, would be restricted to cases of simple mass action kinetics. 
Hence, our work shows that an analytical approach completely distinct from CRNT converges to the 
same strong result, and extends it to a large category of kinetic rate laws. By approaching the question 
of stability with an analytical flavor of SKM, we demonstrate that the dynamics of a well-organized 
metabolic network, even in the presence of nonlinear kinetic rate laws, is not necessarily beyond "pencil 
and paper" methods, and that useful and significant conclusions can be drawn from such analyses. 
Additionally, recent work has shown that a broad category of regulatory and biochemical networks 
can be analyzed using the theory of monotone systems [201 ] . In the future it will be interesting to 
determine whether this theory applies to the cycles discussed here, and conversely how our analytical 
use of SKM may tie into the theory of monotone systems. 

In applying the same SKM techniques to a simple autocatalytic cycle (Figure SJ, we found that 
a combination of analytically derived conditions and computational experiments could highlight the 
dependence of stability on metabolite concentrations, fluxes, and degrees of saturation. In this case, the 
stability of the cycle was not guaranteed, and the analysis suggested possible directions that evolution 
may take in order to steer the network towards an even more robust steady state. While it was possible 
to derive expressions equivalent to Eq. (|5.4p without employing the SKM non-dimensionalization, it 
should be obvious that the algebra required to do so would be complicated and non-trivial assuming 
the rate laws were nonlinear. SKM provided an appealing (but nevertheless rigorous and useful) 
workaround, and still produced the correct characteristic equation for the Jacobian. To date, we are 
not aware of a similar effort to characterize the dynamic properties of autocatalytic cycles for general 
kinetic rate laws. 

Our results suggest that a system does not necessarily require complex forms of regulation in order 
to maintain a stable steady state. In fact, in many cases it seems plausible that the elevated concentra- 
tion of a single metabolite may be sufficient to endow stability. Why is it then that we observe many 
complex forms of allosteric (and genetic) enzyme regulation among biochemical networks? Certainly, 
such regulation may play a crucial role in stabilizing more complex enzymatic networks, but perhaps 
there is more to the story. One possibility may be that these forms of regulation actually control 
the transition of the system from one steady state to another. When conditions in the environment 
change, it may be favorable for a cell to alter its metabolism to react to these changes. Allosteric 
and transcriptional regulation may determine the manifolds along which the transitions between two 
steady states actually occur. 

It is also interesting to consider the possibility that the inherent stability of certain types of path- 
ways may influence the evolution of biochemical networks. One possible link between dynamic stability 
and evolution might be the adaptive fine tuning of kinetic parameters to guarantee stability of a given 
pathway topology. Conversely, we are suggesting that natural selection might gradually amplify the 
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preference for network topologies with an intrinsically high probability of being stable. Such topolo 



gies would be robust and not require substantial fine tuning of kinetic parameters [2 lL |22| . Defining 
the precise parameter ranges and topological structures which promote stability among these pro- 
cesses may enable one to actually search for these common motifs across a broad range of well-studied 
metabolic networks. As more data is collected characterizing the rates of reactions and magnitudes 
of concentrations in metabolism, the possibility of such a search becomes more likely. Searching for 
dynamical behaviors by varying topologies and parameter ranges is quickly becoming a prototypical 
problem in systems biology. Work such as [23l . [24[ demonstrate that solutions to these problems can 
have quite fruitful repercussions in our ability to understand naturally occurring biological systems 
and synthetically design novel genetic and signalling networks. 

Finally, the observation that metabolic cycles have robust steady state properties may have ram- 
ifications in the study of early metabolism. Several studies have suggested that the early stages of 
molecular self-organization and self-reproduction may have involved the emergence of self-sustaining 
proto-metabolic cycles [HQ. Arguments such as those presented in (25[, however, remind us that the 
earliest biological reactions may have taken place in environments which siphoned off the products of 
these reactions. In the case of cycles, these dilutions would have potentially precluded a replenishment 
of substrates in the cycle, and thus prevented a stable steady state. Irrespective of this open question, 
our work suggests that when autocatalytic cycles did first appear, they would have benefitted from 
inherent robustness to environmental perturbations. In general, these cycles may be only a single in- 
stance of a large array of other metabolic topologies which exhibit equally robust stability properties. 
Future studies may address whether these stability properties played a role in establishing the core of 
current cellular metabolism. 
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Appendix A. Structural Kinetic Modeling (SKM) 



The methods used in this paper are based on the SKM framework . SKM is a specific application 
of generalized modeling flo| in which normalized parameters replace conventional kinetic parameters 
(such as V max or Km) in the modeling of metabolic networks. The normalized parameters have a 
direct connection to the original kinetic parameters, but are much easier to work with. As will be 
shown below, these parameters usually have well-defined and extremely limited ranges (e.g. [0,1]), and 
sampling them across this range effectively samples all possible values of the original kinetic parametsrs. 

The goal of SKM is to capture the local stability properties of a biochemical system. In this sense, 



it bridges the gap between genome-scale steady state modeling [ll| and explicit kinetic modeling of a 
metabolic process. To accomplish this, one usually determines the Jacobian J of the system of interest 
and evalutes it at its equilibrium. Assuming knowledge of the form of kinetic rate laws for a CRN, it 
is quite possible to write down the corresponding J. However, in many cases J will be unnecessarily 
complicated and quite difficult to work with. By performing a change of variables, SKM actually 
simplifies the mathematical form of each entry in J. The new entries are in almost all cases easier to 
work with. Much of this material is taken directly from Q and we refer readers to that reference and 
its supplementary materials for more information. 

If we let S be the m-dimensional vector of metabolite concentrations, N be the m x r stoichiometric 
matrix, and v be the r-dimensional vector of reaction rates, then we can describe the dynamics of the 
system with the equation 

§ = N»(S,k) (A.1) 

where w(S,k) denotes that the reaction rates are dependent on both metabolite concentrations S 
and unknown parameters (such as Michaelis-Menten constants) k. If we assume that a non-negative 
steady state S° exists, then we can redefine our system using the definitions 

*■ - M> A- -N-^l a-(x) ~ fA2) 
x i — s o t-'Hf — iV y s o ' — w .(so) ^ ' 

where i — l...m and j = l...r. Now, x is a vector of metabolite concentrations normalized with 
respect to their steady state concentrations and /i represents flux normalized with respect to steady 
state flux values. The matrix A represents the stoichiometric matrix normalized with respect to steady 
state fluxes and steady state metabolite concentrations. We can now rewrite the system of differential 
equations as 

f = A/x(x) (A.3) 

Now, we can write the Jacobian and evaluate it at x° = 1 (which, because of the way x is defined, 
is actually the equilibrium of the system). Calculating the Jacobian, we find that 



= A^M (A4) 
ox 



If we evaluate (|A.4I) at x° = 1 we find 



J x = Ad£ (A.5) 

Note that the equations above were derived without regard to the actual form of the kinetic 
equations that determined the ODE system. The matrix 9£ contains elements which represent the 
degree of saturation of normalized flux fij with respect to normalized substrate concentration Xj. 
In terms of derivatives, each element of 9 represents the degree of change in a flux as a particular 
metabolite is incrementally increased. This is analogous to the concept of elasticity in metabolic 
control analysis [12j . 
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What does the 9 matrix look like? Its columns correspond to each metabolite, and its rows to each 
flux. A non-zero element in the matrix represents the effect a small change in a metabolite j has on 
flux i. In the case of Michaelis-Menten kinetics, this element in the matrix may take values ranging 
from [0,1]. In the case of standard competitive inhibition (e.g. allosteric inhibition by a product), the 
element takes values in [-1,0]. 

To usefully illustrate the meaning of a single 9 element, consider an equation following Michaelis- 
Menten kinetics shown in (|2.4p . First, we write /J-(x), which we recall is the flux normalized by the 
flux at the steady state. Here, So is the concentration of the substrate S at steady state. Then, we 
manipulate the equation by substituting xSo for S, where x is the normalized steady state concentration 
of substrate S. The result is shown in (|2.5j) . 

/ \ Km+s Km + 5*0 7 > 
H\x) — t I g = x— — (A. 7 

«M+f>0 

Finally, we take a derivative with respect to x and evaluate it at x = 1 to obtain 9 in (|2.6p . Notice 
that 9 can only take values between and 1 for any positive value of So ■ 



1 



! , s, 



(A.S 



Km 



It is also possible to approximate the value of 9 for different extremes of saturation. In the regime 
of extremely low substrate concentrations when the enzyme catalyzing the reaction is not saturated 
at all, formation of product can be approximated by 

dP = k 2 E Q S ^ k 2 E S 
dt K M + S~ K m 

which is approximately linear in S. In this case, 9* would be approximately 1. This is identically 
the case for mass-action kinetics, where the rate of product formation is linear with the amount of 
substrate. In contrast, at extremely high substrate concentrations, 

dP k 2 E S ^ 

^-Km~Ts~ 2 ( 0) 

In this case, S has nearly no effect on the rate of product formation, and so 9* is approximately 
0. For all other substrate concentrations, 9 falls between and 1. Note that the above analysis can 
be performed for other forms of kinetic equations including mass-action and Hill-type equations. In 
these cases, the values for 6** are similarly constrained. For example, in the case of Hill kinetics, 9^ 
falls in the range [0,n], where n is the Hill constant. For more information on analyzing other types of 
kinetics, see the supplementary materials in Q. 

The power of SKM comes from the parameterization illustrated above. Each element in 9 has a 
precise correspondence to some combination of kinetic parameters in the original model. However, it 
is far more tractable to study a system using 9 parameters rather than the original kinetic parameters. 
To see this more clearly, recall that in most cases, biochemical kinetic constants are poorly estimated. 
In order to build precise ODE models for biochemical systems, it is usually necessary to actually choose 
values for these constants. While the chosen values may be estimated from experimental measurements, 
hidden dependencies in the model may actually results in non-obvious correlations between parameters 
that can strongly affect the output of the model. 

On the other hand, if one is only interested in the stability of a characterized steady-state, it may 
not be necessary to actually have knowledge of kinetic parameters. Experimental measurements can 
provide data on absolute metabolite concentrations and flux values, and the stoichiometry is often 
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known a priori. Then, one can parameterize the system as shown above and sample many possible 
combinations of 9 parameters. For each unique set of 9 parameters, the stability of the Jacobian is 
determined and recorded. Analysis of the results can lead to several important conclusions such as 
which 9 parameters have the strongest correlation with stability of the entire metabolic system. Thus, 
the benefit of employing SKM over other techniques is that it provides the means to analyze and make 
sense of a large number of possible cases of a metabolic network, rather than a single instance. 

Traditionally, SKM has been used as a computational tool to reveal relationships between the 
stability of a metabolic module and the parameters which govern its kinetics. However, many aspects 
of the SKM parameterization lend themselves to analytical methods. By using 9 parameters, it is no 
longer necessary to decide on the kinetic form (e.g. mass action, Michaelis-Menten, Hill) of the chemical 
reactions taking place. Instead, the problem of studying stability is reduced to finding the signs of the 
roots of a polynomial equation (the characterstic equation of the Jacobian J), whose coefficients are 
determined by the fluxes, absolute metabolite concentrations, and 6's of the system at steady state. 
Moreover, it is in many cases possible to limit the values of 9 (e.g.- if simple Michaelis-Menten kinetics 
are assumed without inhibition, 9 falls between and 1). If the signs of the roots are all negative, 
then the eigenvalues of the Jacobian are negative and the current steady state is stable. 

At first glance, it is unclear whether this will make things easier or harder; values for fluxes and 
concentrations are almost never known. However, our work shows that the absence of this knowledge 
does not make the problem of stability intractable. Instead, we operate with the following idea: we 
assume that the steady state flux vector v and the metabolite concentration vector X are variables, 
and we do not specify them. Instead, we simply determine the stability of the system in terms of these 
variables. Then, we can find trends in stability as the flux and metabolite concentrations change: it 
may be that as a concentration goes up, the stability of the system tends to fall. Isolating these trends 
is the motivation of our analysis. 

It is worthwhile to clarify what it means to change a steady-state parameter like a flux or a 9. At 
first glance, it may not be obvious that it is even possible to alter such a parameter. In other words, 
suppose that at steady state, the concentration of metabolite A in some network is 1. What does 
it mean to increase the concentration of Al If we do so without altering any other parameters, the 
metabolic network will in all likelihood not be in steady state any longer. However, there may exist 
(and in fact for the systems described in this paper there does exist) another set of concentrations, 
fluxes, and 9's which does lead to a steady state at the increased concentration of A. The problem, 
again, is that we do not know their actual values. However, we may be able to infer, based on changes 
in the coefficients of the characteristic equation, that if wc were to sample a large number of unique 
instances of the Jacobian for both concentrations of A, the higher concentration may actually be more 
likely to be stable. This is in fact the method we will use to study autocatalytic cycles later in the 
paper. However, in the next two sections, we will demonstrate that for certain types of cycles, it is 
not necessary to go to such great lengths in order to characterize stability; it will come much more 
naturally. 

Appendix B. Linear Chains 

Up to this point, linear metabolic networks have been ignored. The reason for this is that, using 
the same analysis as above, it is trivial to show such a network is stable. In this appendix, we provide 
a short example illustrating this point. Consider the metabolic network shown in Figure 3 with flux 
F flowing through it and observed steady states (A°B°). 

If we follow the same analysis as before (writing A and 9 and calculating the Jacobian J), we will 
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Figure B.6: A linear chain 

The eigenvalues of the Jacobian are then and both of which are negative. Therefore, the 
system is stable. 

Appendix C. Cofactors 

The reaction system shown below is the complete, more general version of the simple system 
analyzed above. 

Now, we can write a system very similar to the one in Section 2. Consider the example shown in 
the figure below. 
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Figure C.7: A complete metabolic cycle with cofactors 

The generalized form of A and 6 is shown below. Note that we make no further assumptions about 
the steady states of the system: A" = (X^X®, ■ ■ ■ ,X^) and the flux through the system remains 
equal to F. With N metabolites and N + 3 reactions, A is N x (N + 3) and 6 is (N + 3) x N. 
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Note that because the cofactors come in a conserved pair, one of them can be removed from the 
system. This modifies the entry in the bottom right element of the Jacobian. 
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Note that the last row of A is omitted because the cofactors come as a conserved pair. Furthermore, note that the bottom right element of 
replaced with a negative element in order to account for this conservation. Aftering factor out F, the Jacobian becomes 
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Next, we evaluate take the determinant | J — A/| and take the Laplace expansion over the last column. This yields three smaller determinants. 
The first is ~^o +2 — ^rr 2 multiplied by the determinant of: 
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The second matrix resulting from the Laplace expansion of the original Jacobian is (— 1) 
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If another Laplace expansion is now taken along the N th row, then the resulting matrix is 
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The resulting term for this determinant is 
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The third matrix resulting from the Laplace expansion is ( — 1) W +- R + 2 ^"+ 2 multiplied by the determinant of 
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Taking another Laplace expansion over the N th column, then the resulting matrix is 
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The resulting term for this determinant is 



Now we can sum all terms to get 
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The terms on the bottom two lines will cancel with terms from the top two lines, resulting in a polynomial equation of order N+ 1 containing 



nonzero coefficients with exactly the same sign. Furthermore, a constant A term will remain, indicating that there is no zero eigenvalue. 



